
#Replication code for HGZ, ZC & RS (2025)
install.packages("mice")
install.packages("lavaan")
install.packages("car")
install.packages("apaTables") 
library(apaTables)
library(mice)
library(foreign)
library(dplyr)
library(lavaan)
library(lm.beta)
library(car)
############################################
#### Wave 1 - cross-sectional model ########
############################################
df <- read.spss ("/Users/clairecheng/Desktop/Repression/HSSC_submission/RR1/replication_data/HSSC_replication data/w1_HSSC_rep.sav", use.value.labels = TRUE, to.data.frame=TRUE)

####Cross-sectional Model 
mydf<-df %>% dplyr::select(age_r,sex_r, Q3, Q4, RACE_r, Q368r, Q368, ideo, tra, smn, opd, obst,cb)
mydf$age_r<-as.numeric(mydf$age_r)
mydf$Q3<-as.numeric(mydf$Q3)
mydf$Q4<-as.numeric(mydf$Q4)
mydf$sex_r<-as.numeric(mydf$sex_r)
mydf$RACE_r<-as.numeric(mydf$RACE_r)
mydf$Q368<-as.numeric(mydf$Q368)

#impute
set.seed(432)
mydf_mice1 <- mice(mydf,m=1)
df.imp<- mice::complete(mydf_mice1)

####Table 1 (CM as DV)
block3 <- lm(cb~age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+obst, data=df.imp)
summary(lm.beta(block3),standardized=TRUE)
apa.reg.table(block3, filename="apa_regression_table.doc")


####Table 2 (PI as DV)
mydf<-df %>% dplyr::select(age_r,sex_r, Q3, Q4, RACE_r, Q368r, Q368, ideo, tra, smn, opd, obst,cb,pr)
mydf$age_r<-as.numeric(mydf$age_r)
mydf$Q3<-as.numeric(mydf$Q3)
mydf$Q4<-as.numeric(mydf$Q4)
mydf$sex_r<-as.numeric(mydf$sex_r)
mydf$RACE_r<-as.numeric(mydf$RACE_r)
mydf$Q368<-as.numeric(mydf$Q368)

#impute
set.seed(432)
mydf_mice1 <- mice(mydf,m=1)
df.imp<- mice::complete(mydf_mice1)
summary(mydf_mice1)
summary(df.imp)

block3 <- lm(pr~age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+obst+cb, data=df.imp)
summary(lm.beta(block3), standardized=TRUE)
apa.reg.table(block3, filename="apa_regression_table.doc")

####Table 3 (Mediations)
set.seed(1234)
sem_model = '
  cb ~ a*obst +age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd
  pr ~ c*obst + age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd + b*cb
 
  # direct effect
  direct := c
 
  # indirect effect
  indirect := a*b
 
  # total effect
  total := c + (a*b)
'

model_sem = sem(sem_model, data=df.imp, se='boot', bootstrap=1000)
fit<-summary(model_sem, standardized = TRUE, ci=TRUE)


################################
####Wave 2 - Lagged model ######
################################
df <- read.spss ("/Users/clairecheng/Desktop/Repression/HSSC_submission/RR1/replication_data/HSSC_replication data/w2_HSSC_rep.sav", use.value.labels = TRUE, to.data.frame=TRUE)

mydf<-df %>% dplyr::select(age_r,sex_r, Q3, Q4, RACE_r, Q368r, Q368, ideo, tra, smn, opd, obst,cb2)
mydf$age_r<-as.numeric(mydf$age_r)
mydf$Q3<-as.numeric(mydf$Q3)
mydf$Q4<-as.numeric(mydf$Q4)
mydf$sex_r<-as.numeric(mydf$sex_r)
mydf$RACE_r<-as.numeric(mydf$RACE_r)
mydf$Q368<-as.numeric(mydf$Q368)

#impute
set.seed(432)
mydf_mice1 <- mice(mydf,m=1)
df.imp<- mice::complete(mydf_mice1)

####Table 1 (CM as DV)
block3 <- lm(cb2~age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+obst, data=df.imp)
summary(lm.beta(block3),standardized=TRUE)
apa.reg.table(block3, filename="apa_regression_table.doc")

####Table 2 (PI as DV)
mydf<-df %>% dplyr::select(age_r,sex_r, Q3, Q4, RACE_r, Q368r, Q368, ideo, tra, smn, opd, obst,cb,pr2)

mydf$age_r<-as.numeric(mydf$age_r)
mydf$Q3<-as.numeric(mydf$Q3)
mydf$Q4<-as.numeric(mydf$Q4)
mydf$sex_r<-as.numeric(mydf$sex_r)
mydf$RACE_r<-as.numeric(mydf$RACE_r)
mydf$Q368<-as.numeric(mydf$Q368)

#impute
set.seed(432)
mydf_mice1 <- mice(mydf,m=1)
df.imp<- mice::complete(mydf_mice1)

block3 <- lm(pr2~age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+obst+cb, data=df.imp)
summary(lm.beta(block3), standardized=TRUE)
apa.reg.table(block3, filename="apa_regression_table.doc")

####Table 3 (Mediations)
set.seed(1234)
sem_model = '
  cb ~ a*obst +age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd
  pr2 ~ c*obst + age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd + b*cb
 
  # direct effect
  direct := c
 
  # indirect effect
  indirect := a*b
 
  # total effect
  total := c + (a*b)
'
model_sem = sem(sem_model, data=df.imp, se='boot', bootstrap=1000)
fit<-summary(model_sem, standardized = TRUE, ci=TRUE)  


#######################################
#### Wave 2 - Autoregressive model ####
######################################

df <- read.spss ("/Users/clairecheng/Desktop/Repression/HSSC_submission/RR1/replication_data/HSSC_replication data/w2_HSSC_rep.sav", use.value.labels = TRUE, to.data.frame=TRUE)

####Table 1 (CM as DV)
mydf<-df %>% dplyr::select(age_r,sex_r, Q3, Q4, RACE_r, Q368r, Q368, ideo, tra, smn, opd, obst,cb,cb2)
mydf$age_r<-as.numeric(mydf$age_r)
mydf$Q3<-as.numeric(mydf$Q3)
mydf$Q4<-as.numeric(mydf$Q4)
mydf$sex_r<-as.numeric(mydf$sex_r)
mydf$RACE_r<-as.numeric(mydf$RACE_r)
mydf$Q368<-as.numeric(mydf$Q368)

#impute
set.seed(432)
mydf_mice1 <- mice(mydf,m=1)
df.imp<- mice::complete(mydf_mice1)

block3 <- lm(cb2~cb+age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+obst, data=df.imp)
summary(lm.beta(block3),standardized=TRUE)
apa.reg.table(block3, filename="apa_regression_table.doc")

####Table 2 (PI as DV)
mydf<-df %>% dplyr::select(age_r,sex_r, Q3, Q4, RACE_r, Q368r, Q368, ideo, tra, smn, opd, obst,cb,pr2,pr)

mydf$age_r<-as.numeric(mydf$age_r)
mydf$Q3<-as.numeric(mydf$Q3)
mydf$Q4<-as.numeric(mydf$Q4)
mydf$sex_r<-as.numeric(mydf$sex_r)
mydf$RACE_r<-as.numeric(mydf$RACE_r)
mydf$Q368<-as.numeric(mydf$Q368)


#impute
set.seed(432)
mydf_mice1 <- mice(mydf,m=1)
df.imp<- mice::complete(mydf_mice1)

block3 <- lm(pr2~pr+age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+obst+cb, data=df.imp)
summary(lm.beta(block3), standardized=TRUE)
apa.reg.table(block3, filename="apa_regression_table.doc")

####Table 3 (Mediations)
set.seed(1234)
sem_model = '
  cb ~ a*obst +age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+pr
  pr2 ~ c*obst + age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+pr + b*cb
 
  # direct effect
  direct := c
 
  # indirect effect
  indirect := a*b
 
  # total effect
  total := c + (a*b)
'
model_sem = sem(sem_model, data=df.imp, se='boot', bootstrap=1000)
fit<-summary(model_sem, standardized = TRUE, ci=TRUE)  



####Sensitivity analysis
####sensitivity analysis on lagged 
install.packages("mediation")
library(mediation)
b <- lm(cb ~ obst + age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd, data=df.imp)
c <- lm(pr2 ~ obst + cb + age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd, data=df.imp)
med.cont <- mediate(b, c, treat="obst", mediator="cb", sims=100)

sens.cont <- medsens(med.cont, rho.by=.05)
summary(sens.cont)

par.orig <- par(mfrow = c(2,2))
plot(sens.cont, main="Lagged", ylim=c(-.2,.2))
plot(sens.cont, sens.par="R2", r.type="total", sign.prod="positive")
plot(sens.cont, sens.par="R2", r.type="total", sign.prod="negative")#when the unobserved confounder affect M and DV in opposite directions. 
par(par.orig)

####sensitivity analysis on autogressive 
#select the imputed df in the above block  
b1 <- lm(cb ~ obst + age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+pr, data=df.imp)
c1 <- lm(pr2 ~ obst + cb + age_r+sex_r+Q3+Q4+RACE_r+Q368r+Q368+ideo+tra+smn+opd+pr, data=df.imp)
med.cont <- mediate(b1, c1, treat="obst", mediator="cb", sims=100)

sens.cont <- medsens(med.cont, rho.by=.05)
summary(sens.cont)

par.orig <- par(mfrow = c(2,2))
plot(sens.cont, main="Autoregressive", ylim=c(-.2,.2))
plot(sens.cont, sens.par="R2", r.type="total", sign.prod="positive")
plot(sens.cont, sens.par="R2", r.type="total", sign.prod="negative")#when the unobserved confounder affect M and DV in opposite directions. 
par(par.orig)




